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A system of coupled oscillators can exhibit a rich variety of dynamical behaviors. When we investigate the 
dynamical properties of the system, we first analyze individual oscillators and the microscopic interactions 
between them. However, the structure of a coupled oscillator system is often hierarchical, so that the 
collective behaviors of the system cannot be fully clarified by simply analyzing each element of the system. 
For example, we found that two weakly interacting groups of coupled oscillators can exhibit anti-phase 
collective synchronization between the groups even though all microscopic interactions are in-phase 
coupling. This counter-intuitive phenomenon can occur even when the number of oscillators belonging to 
each group is only two, that is, when the total number of oscillators is only four. In this paper, we clarify the 
mechanism underlying this counter-intuitive phenomenon for two weakly interacting groups of two 
oscillators with global sinusoidal coupling. 

A system of coupled oscillators provides abundant examples of dynamical behaviors including synchron- 
ization phenomena' Among them, collective synchronization emerging from coupled phase oscillators 
has been widely investigated not only for globally coupled systems but also for complex network sys- 
tems'^"". Furthermore, the dynamical behaviors exhibited by interacting groups of globally coupled phase 
oscillators have been intensively investigated"*"^'. The appearance of the Ott-Antonsen ansatz^**"™ has consid- 
erably facilitated theoretical investigations on interacting groups of noiseless nonidentical phase oscillators with 
global sinusoidal coupling. In addition, interacting groups of globally coupled phase oscillators as well as a system 
of globally coupled phase oscillators have been experimentally realized using electrochemical oscillators^''^^, 
discrete chemical oscillators'''''''', and mechanical oscillators''^""'. 

To study the phase synchronization between macroscopic rhythms, we recently formulated a theory for the 
collective phase description of macroscopic rhythms emerging from coupled phase oscillators for the following 
three representative cases: (A) phase coherent states in globally coupled noisy identical oscillators''^ (B) 
partially phase-locked states in globally coupled noiseless nonidentical oscillators''", and (C) fully phase-locked 
states in networks of coupled noiseless nonidentical oscUlators"". The theory enables us to describe the dynamics 
of a macroscopic rhythm by a single degree of freedom called the collective phase. Accordingly, different 
mathematical treatments were required for the physical situation in each case. The keystone of the collective 
phase description method for each case is the following: (A) the nonlinear Fokker- Planck equation^, (B) the Ott- 
Antonsen ansatz^**"'", and (C) the Laplacian matrix'*"". Here, we note that there exist several investigations*^"'"' 
related to case (C). 

In Ref 39 for case (A) and Ref 40 for case (B), we investigated the phase synchronization between collective 
rhythms of globally coupled oscillator groups. In particular, the collective phase coupling function, which 
determines the dynamics of the collective phase difference between the groups, was systematically analyzed 
for sinusoidal coupling functions. As a result, for both cases, we found counter-intuitive phenomena in which 
the groups can exhibit anti-phase collective synchronization in spite of microscopic in-phase external coupling 
and vice versa. 

In this paper, using the collective phase description method developed in Ref. 41, we study the phase syn- 
chronization between collective rhythms of coupled oscillator groups for case (C). We analytically derive the 
collective phase coupling function for two weakly interacting groups of two oscUlators with global sinusoidal 
coupling (see Fig. 1). We thereby demonstrate counter- intuitive phenomena simOar to those found in cases (A) 
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Figure 1 | Schematic diagram of two weakly interacting groups of two 
oscillators with global coupling. The microscopic internal and external 
couplings are represented by the solid and dotted arrows, respectively, 
whereas the self-coupling is not shown. The phase of the j-th oscillator in 
the (7-th group is denoted by i/i . 



and (B), that is, effective anti-phase (in-phase) collective synchron- 
ization with microscopic in-phase (anti-phase) external coupling. 
Therefore, this paper and Refs. 39, 40 are mutually complementary 
and together provide a deeper understanding of the collective phase 
synchronization phenomena. 

Results 

This section is organized as follows. First, we formulate the collective 
phase description of fully locked states with an emphasis on the 
collective phase coupling function. Second, we analyze weakly inter- 
acting groups of globally coupled two phase oscillators. Third, we 
perform further analytical calculations for the case of sinusoidal 
phase coupling. Fourth, we illustrate the collective phase coupling 
function for several representative cases. Fifth, we demonstrate col- 
lective phase synchronization by direct numerical simulations. 
Finally, we consider interacting groups of weakly coupled Stuart- 
Landau oscillators. 

Collective phase description of fully locked states. We consider 
weakly interacting groups of coupled noiseless nonidentical phase 
oscillators described by the following equation: 



(1) 



i:=l 



for; = 1 , • • • , N and ((7, t) = (1, 2), (2, 1), where f"' {t)eS^ is the phase 
of the j-th oscillator at time t in the ff-th group consisting of N oscillators 
and cOj is the natural fi^equency of the j-th phase oscillator. The second 
term on the right-hand side represents the microscopic internal 
coupling within the same group, while the third term represents the 
microscopic external coupling between the different groups. The 
characteristic intensity of the external coupling is given by 6>0. 
When the external coupUng is absent, i.e., £ = 0, Eq. (1) is assumed to 
have a stable fiiUy phase-locked collective oscillation solution'"'''"'™ 



where (t)eS^ is the collective phase at time t for the cr-th group, Q 
is the collective frequency, and the constants ij/j represent the relative 
phases of the individual oscillators for the fiiUy phase-locked state. 



When the external coupling is sufficiently weak, i.e., e<d, each 
group of oscillators obeying Eq. (1) is always in the near vicinity of 
the fuUy phase-locked solution (2). Therefore, we can approximately 
derive a collective phase equation in the following form'*': 



where the collective phase coupling function is given by 

(©(-) _ 0(^)) = j2fl f®'"' - + ^; - 

j=l k=l 



(3) 



(4) 



Here, L^* is the left zero eigenvector of the Jacobi matrix Lji^ at the 
fully phase-locked collective oscillation solution defined in Eq. (2). 
The Jacobi matrix Lji^ is given by 



(l-<5;0r;4^,-^J, (5) 



which is a Laplacian matrix'''". That is, the Jacobi matrix L^j- pos- 
sesses the following property for each j: 2_^k- 1 ^i*: ^ ^- ^1- 
have used the Kronecker delta 5jt and derivative notation 
= dYjk{(j>) /d<j>. We also note that the Jacobi matrixiy^^ defined 



in Eq. (5) is generally asymmetric and weighted. Using the (j, j)- 
cofactor of the Jacobi matrix and the summation over the index j, i.e.. 



Mj= det L{j,j), 



M-- 



(6) 



the left zero eigenvector U* of the Jacobi matrix that takes the form of 
the Laplacian matrix can be generally written in the following form'""*': 



. M,- 



(7) 



In Eq. (6), the matrix Lij,}) is the Jacobi matrix L with the j-th row and 
column removed, and the cofactor Mj is equal to the sum of the weights 
of all directed spanning trees rooted at the node j according to the 
matrix tree theorem^''^^. Finally, we note that the collective phase 0*°' 
can be written in the following form'": 



(8) 



under the linear approximation of the isochron'"^''*'". 



Interacting groups of globally coupled two phase oscillators. We 

here analyze globally-coupled two-osciUator systems using the 
collective phase description method for fully locked states. We first 
consider weakly interacting groups of globally coupled phase 
oscillators. That is, the microscopic internal and external coupling 
functions are given by 



I;* 



(9) 



n{f-<l^) = ],^''{^f'-^)- (10) 

(2) In this global coupling case, Eq. (1) is written in the following form: 



(11) 
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We further focus on the case in which the number of oscillators 
within each group is two, i.e., N = 2; a schematic diagram of the 
case is shown in Fig. 1. In this case, the internal dynamics for each 
group, i.e., Eq. (11) with e = 0, is described as follows: 



^j(f) = aji + ir(o) + ir(^i-^2), 



(12) 



(13) 



where we dropped the group index a for simplicity. From Eqs. (12) 
and (13), we obtain the following equation by subtraction: 

|A^(O=Affl+ir(A0)-ir(-A^), (14) 

where the phase difference A^(f) and frequency mismatch Aw are 
defined as 



A(j>{t)=(t>i{t)—(f>2{t), Aaj = oJi— C02- 



(15) 



Now, we assume that Eq. (14) has a fuUy phase-locked collective 
oscillation solution. The phase difference of the stable phase- 
locked solution, Aij/ = ipi — i/'2, is determined by the following 
equation: 



A«+ir(A^)-ir(-A^) = o. 



(16) 



Using the phase difference Ai// obtained from Eq. (16), the collective 
frequency Q is written in the following form: 



n=coi + ^r(o) + ^r(A,A) 
= «2+^r(o) + ^r(-A^). 



(17) 



(18) 



(19) 



For these globally-coupled two-oscillator systems, the Jacobi matrfac 
L defined in Eq. (5) is given by 

^^1/ r'iAik) -r'(A,A)\ 
2V-r'(-A^) r'(-A,A)/ 

Therefore, the cofactors of the Jacobi matrix are given by 

r'(-AiA) r'(Ai//) 

Mi = — M2= ^ , 

2 2 

and the sum of these cofactors is written as 

M = Mi+M2=— ^ (20) 

As found from Eq. (7), using these cofactors and the sum, Eq. (19) 
and Eq. (20), the left zero eigenvector U* is obtained as 

Ml r'(-AiA) 



^1= 



M r'(-A(/^)-Hr'(AiA)' 
Mj r'(AiA) 



(21) 



^2 M r'(-AiA) + r'(Ai//)' 

Finally, we note that the Jacobi matrix L possesses not only the zero 
eigenvalue but also the following non-zero eigenvalue: 

r'(-A^) + r'(A^) 



X=M = 



22 



When the external coupling intensity is sufficiently small compared 
to the absolute value of this non-zero eigenvalue, i.e., e^|A|, the 
collective phase description is valid"". 



Analytical formulas for the case of sinusoidal phase coupling. We 

here consider the case of sinusoidal phase coupling functions for 
both microscopic internal and external couplings. First, the 
microscopic internal phase coupling function is given by 

VITjk{^) = T{^) = -sin{^ + ^.), |0!|<|, (23) 

which is in-phase coupling (i.e., attractive). By substituting Eq. (23) 
into Eq. (16), the phase difference of the fuUy phase-locked state is 
obtained as 



sin(Ai//) = !/, 11 E 



Act) 



I<1, 



(24) 



which indicates that the fuUy phase-locked solution emerge from a 
saddle-node bifurcation and exists under the condition of | At;j| < cos 
a. Owing to the in-phase coupling, i.e., Eq. (23), one solution of | Ai//] 
< nl2 is stable, and the other solution of > nil is unstable. 
Hereafter, the fully phase-locked solution indicates the stable one, 
jAi/'l < nil. Substituting Eqs. (23) and (24) into Eq. (17), we obtain 
the collective frequency Q. as 



i2 = t;ji — i^cos a- 



1 + 



(25) 



Similarly, substituting Eqs. (23) and (24) into Eq. (19), we obtain the 
cofactors as follows: 



Mi = 



cos d — r] sm a 



(26) 



- cos a + )7 sm a. 



which yield M= — — cos a. From Eqs. (21) and (26), the left 
zero eigenvector U* is thus written as 



1 / n tan a 



(jtana 



(27) 



In addition, the non-zero eigenvalue 1 defined in Eq. (22) is obtained 
as 



X=- \J\-n^ cos a. (28) 
Next, the microscopic external phase coupling function is given by 



\ = - sm{(l) + p), 



(29) 



which can be either in-phase coupling (i.e., attractive) under the 
condition of |/?| < nil or anti-phase coupling (i.e., repulsive) 
under the condition of > nil. By plugging Eqs. (24), (27), and 
(29) into Eq. (4), the collective phase coupling function takes the 
following form: 



P''^(0) = -psin(©-h^), 



(30) 



where the complex number with modulus p and argument 3 is given 
by 



pe"' = 



COS p- 
sin 



,^2 



:tan a. sm 



(31) 



r 



:tan a cos p 



This formula is the main result of the present paper. It determines the 
collective phase coupling function for two weakly interacting groups 
of two oscillators with global sinusoidal coupling. The coupling type 
can be found from the real part, i.e., 
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p COS (5 = 



1 



2 



ztan a sin /J 



(32) 



where p cos (5 > 0 and p cos (5 < 0 indicate in-phase and anti-phase 
coupHngs, respectively. Finally, we note that Eq. (32) possesses origin 
symmetry in the a-jS plane. 

Type of the collective phase coupling function for representative 
cases. We here study the type of the collective phase coupling 
function for the following five representative cases. 



(i) 



(ii) 



(iii) 



(iv) 



The first case is = 0, which indicates that two oscillators 
within each group are identical, i.e., Acu = 0. Substituting 
(7 = 0 into Eq. (31), we obtain the following result: 



(33) 



That is, the collective phase coupling function is the same as 
the microscopic external phase coupling function, i.e., P" (©) 
= r" (0) = -sin(0 + P). 

The second case is — 1, which indicates the proximity of 
the saddle-node bifurcation point, i.e., the onset of fully 
phase-locked collective oscillation. Substituting — 1 into 
Eq. (31), we obtain the following result: 



pe 



: tan a sin /? 



: tan a cos P 



(34) 



For the case of | | ^ 1 (excluding (X = 0), the amplitude of the 
collective phase coupling becomes infinity, i.e., p Here, 
we note that this property for the fully phase-locked states is 
quite different from those for phase coherent states and par- 
tially phase-locked states^'^''*". For the latter two states, the 
amplitude of the collective phase coupling is finite at the onset 
of collective oscillations. This difference in the properties 
results from the difference of bifurcations. The fuUy phase- 
locked states emerge from saddle-node bifurcations as men- 
tioned above, whereas the phase coherent states and partially 
phase-locked states emerge from supercritical Hopf bifurca- 
tions'''"". 

The third case is a = 0, which yields a microscopic antisym- 
metric internal coupling function. For this case, r\ = Acu. 
Substituting a = 0 into Eq. (31), we obtain the following 
result: 



(35) 



That is, the phase shift 5 of the collective phase coupling 
function is the same as the phase shift [i of the microscopic 
external phase coupling function. 

The fourth cases are special values of p. Substituting P = 0, 
±7!:, ±n/2 into Eq. (31), we obtain the following results: 



pe 



:tana 



(36) 



P=±K, 



pe-' 



K 

2 ' 



pe 



+ - 

i 
2 



ztana 



(37) 



= tano! 



(38) 



(v) 



For microscopic antisymmetric external coupling functions, 
i.e., p = 0, ±n, the type of the collective phase coupling 
function coincides with that of the microscopic external 
coupling function. In contrast, for microscopic symmetric 
external coupling functions, i.e., p = ±n/2, the type of the 
collective phase coupling function is determined by the sign 
of the microscopic internal coupling parameter a. 
The fifth case is P = a, which indicates that the microscopic 
external coupling has the same phase shift as the microscopic 
internal one. Substituting P = a into Eq. (31), we obtain the 
following result: 



/i = a, pe" 



1 



2 cos a 



(39) 



From the condition of |a| < 7c/2, both microscopic internal 
and external coupling functions are in-phase coupling. 
However, the type of the collective phase coupling function 
is anti-phase coupling under the following condition: 



2 n 

cos a < 



(40) 



For the case of jijl ^ 1, the above condition becomes cos^ a. < 
1, which is satisfied for all a except for a = 0. 

Collective phase synchronization between two interacting groups. 

Now, we study counter-intuitive cases under the condition of = 
3/4. The type of the collective phase coupling function is shown in 
Fig. 2, where the solid curves are determined by Eq. (32), i.e., p cos 8 
= 0. Here, we note that the type of the collective phase coupling 
function can be different from that of the microscopic external 
phase coupling function. Two sets of parameters, which were used 
in Fig. 3, are also shown in Fig. 2. 

Two groups of two-oscUlators exhibiting phase-locked states were 
separately prepared with their corresponding phases being nearly 
identical. Then, these states were used as the initial condition in 
Fig. 3(a). In spite of the microscopic in-phase external coupling, 

P = 3n/8, the external phase difference ji^'/' — i^i^' j approached n 

after some time; this indicates anti-phase collective synchronization 
between the groups. In contrast. Fig. 3(b) shows in-phase collective 
synchronization between the groups in spite of the microscopic anti- 
phase external coupling, p = —5n/&. 

Interacting groups of weakly coupled Stuart-Landau oscillators. 

We further consider interacting groups of globally coupled Stuart- 
Landau oscillators described by the following equation: 
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-nl2 0 kI2 



[T) = 3/4] a 

Figure 2 | Effective type of phase coupling between collective rhythms of 
fully locked oscillator groups with or 6 ( — Jil2,nl2),fi& [— ;r, ^],and^ = 
3/4. The solid curves are analytically determined by Eq. (32), i.e., p cos <5 = 
0. The filled circle (•) indicates y. = ft = 3n/8 corresponding to Fig. 3(a) 
and Fig. 4. The times sign (X) indicates a = 3nl8 and ft = —5nl8 
corresponding to Fig. 3(b). 



Wj"^ it) = {l + ibj) Wj"^ -{l + ic2)\ Wj"^ I Wj"'' 

k=i k=i 

forj= 1, ■ ■ ■ ,N and (p, z) = (1, 2), (2, 1), where wj''\t)eC^ is the 
complex amplitude of the j-th limit-cycle oscillator at time t in the a- 
th group consisting of N oscillators. The first and second terms on the 
right-hand side represent the intrinsic dynamics of each oscillator, 
the third term represents the microscopic internal coupling within 
the same group, and the fourth term represents the microscopic 
external coupling between the different groups. When the internal 
and external couplings are sufficiently weak compared to the 
absolute value of the amplitude Floquet exponent, we can 
approximately derive a phase equation in the following form^: 

/:— 1 



where the parameters of phase oscillators are given by 

a}j = bj-C2, (43) 



PKe-'=Kil + ic2){l-ici), (44) 



P,e'"=}{l + ic2)il-ic3). (45) 

The phase of each Stuart- Landau oscillator is given by the following 
equation' (f) = arg W — C2 In | W|. Here, we focus on the case in 
which the number of oscillators within each group is two, i.e., N =2. 
Using the following constants, r = 0.01 and a = 3n/8, the parameters 
of the Stuart-Landau oscillators are fixed zXK = ] = r cos(a), Ci = C3 
= 0, C2 = tan(a), bi = C2 + 3r cos(a)/4, and ^2 = Ci- Under these 
conditions, the parameters of the p base oscillators are obtained as 
= Pj = r = 0.01, a = P = a = 3n/8, coi = 3r cos(a)/4, and 0J2 = 0, 
which correspond to the parameters in Fig. 3(a). In particular, we 
note that )/ = {Aco)/{PgCos a) = (3r cos(a)/4)/(r cos(a)) = 3/4. The 
external coupling intensity is fixed at £ = 0.001. The direct numerical 
simulation result of Eq. (41) is shown in Fig. 4. Similarly to Fig. 3(a), 
Fig. 4 shows anti-phase collective synchronization between the 
groups in spite of the microscopic in-phase external coupling. 

Discussion 

In this paper, we considered the phase synchronization between 
collective rhythms of fully locked oscillator groups, clarified the rela- 
tion between the collective phase coupling and microscopic external 
phase coupling functions, analytically determined the type of the 
collective phase coupling function for weakly interacting groups of 
two oscillators with global sinusoidal coupling, and demonstrated 
that the groups can exhibit anti-phase (in-phase) collective syn- 
chronization in spite of microscopic in-phase (anti-phase) external 
coupling. The theoretical predictions were successfully confirmed by 
direct numerical simulations of the phase oscillator model and 
Stuart-Landau oscillator model. 

In Refs. 39, 40, we investigated the phase synchronization between 
collective rhythms of globally coupled oscillator groups under two 
typical situations: phase coherent states in the noisy identical case'' 
and partially phase-locked states in the noiseless nonidentical case''". 
In particular, we found the counter- intuitive phenomena similar to 
the results in this paper. That is, weakly interacting groups can 
exhibit anti-phase collective synchronization in spite of microscopic 
in-phase external coupling and vice versa. Here, we note that these 




0.0 1.0 2.0 3.0 0.0 1.0 2.0 3.0 

[/3= 37r/81 r/lO" [i3 = -5;r/81 f/ lO" 



Figure 3| Interacting groups of phase oscillators (see Methods). Time evolution of the internal and external phase differences, i.e., —(l>2^ | and 
jijdj —ij>[ ' I . The other internal and external phase differences are approximated as jijd'j^' ~<i^r'| — | ji^j'' ""^2^' | — —ijif |, respectively. 

The collective phase difference is approximated as the external phase difference, i.e., j©''' —0'^' | ~ ji/ij'' — i/ij^' |. The parameters are a = in/S, coi = 3 
cos((x)/4, CO2 = 0, and e = 0.001. (a) Effective anti-phase collective synchronization with microscopic in-phase external coupling, ft = 3n/8. (b) Effective 
in-phase collective synchronization with microscopic anti-phase external coupling, ft = — Stt/S. 
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0.0 1.0 2.0 3.0 -1.0 0.0 1.0 

[SL] f/lO" ReW 



Figure 4 | Interacting groups of weakly coupled Stuart-Landau (SL) oscillators (see Methods). Effective anti-phase collective synchronization with 
microscopic in-phase external coupling. The parameters areK= /= rcos(fl), Ci = = 0, C2 = tan(fl), bi = C2 + 3rcos(fl)/4, = C2, andf = 0.001, where r 
= 0.01 and a = 3k/S. (a) Time evolution of the internal and external phase differences, i.e., jijfij'' —(l>2^ | and jijdj'' — (l>i |. (h) Snapshot of the asymptotic 
state of individual oscillators, i.e., w{'', Wj'', wj^', and w'^K 



three papers considered different physical situations and utilized 
different mathematical methods, but arrived at the similar counter- 
intuitive phenomena. 

We also remark that fully phase-locked states emerge from a finite 
number of oscillators'"'; even two is possible as actually studied in 
this paper. In contrast, phase coherent states and partially phase- 
locked states emerge from a large population of oscillators^; the 
number of oscUlators is infinite in theory. From this point of view, 
fully phase-locked states can be more easily realized in experiments 
such as electrochemical oscillators'" '^, discrete chemical oscilla- 
tors'^'''*, and mechanical oscillators'^"'. We hope that the counter- 
intuitive phenomena studied in this paper, i.e., effective anti-phase 
(in-phase) collective synchronization with microscopic in-phase 
(anti-phase) external coupling, wall be experimentally confirmed in 
the near future and that the formula (31) wiU help in such 
experiments. 

Finally, we emphasize that collective synchronization between 
interacting groups of coupled oscillators cannot be fully clarified 
by simply analyzing microscopic interactions between individual 
oscillators. In particular, microscopic in-phase (anti-phase) external 
coupling does not necessarily lead to in-phase (anti-phase) collective 
synchronization. As clarified in this paper, counter-intuitive phe- 
nomena can occur even when the number of oscillators belonging 
to each group is only two. We hope that the analytical results for the 
simple cases studied in this paper wUl provide an insight into more 
complex cases. 

Methods 

Numerical method for Fig. 3. We applied an explicit Euler scheme with a time step 
A( = 0.01 forEq. (11) withEqs. (23) and (29). The parameters are JV= 2, a = 37t/8, coi 
= 3 cos(o()/4, CO2 = 0, and c = 0.001 with (a) [i = in/& or (b) /? = -57t/8. The initial 
values are (j>\^\o) = arcsin(3/4) and ,^<''(0) =0 with (a) i^j^'(O) ~(^]''(0) or (b) 
,^f'(0)^^("(0)-H7r. 

Here, we note the accuracy and stability of the numerical method. On the right- 
hand side of Eq. (11), the first and second terms represent the internal dynamics of 
0{ 1 ) whUe the third term represents the external coupling of 0( c ) . When the external 
coupling is sufficiently weak, the smallest time scale in Eq. (11) is 0(1). Therefore, the 
explicit Euler scheme with the time step At — 0.01, which we used for the sake of 
simplicity and efficiency, is sufficiently accurate and stable under the parameter 
condition. 

Numerical method for Fig. 4. We applied an explicit Euler scheme with a time step 
At = 0.01 for Eq. (41). The phase of they'-th oscillator at time t in the c-th group was 

obtained by ^j''^ (f) - argWj''*(f)-C2 ln|wj'''(f)|. The parameters areN- 2,iC--7 

= r cos((i), Ci = = 0, C2 = tan(a), bi = C2 + 3r cos((3)/4, /?2 ~ '^2> snd c — 0.001, 
where r = 0.01 and a = 3n/S. The initial condition is given by 

wj^'^O) - exp(/(j^j^*(0)), where <f)^^\o) - arcsin(3/4), ^^^'(0) -0, and 
•^f (O)-^f' (0). 



We also note the accuracy and stability of the numerical method. On the right- 
hand side of Eq. (41), the first and second terms represent the oscillator dynamics of 
0(1), the third term represents the internal coupling of 0{K), and the fourth term 
represents the external coupling of 0{cj). When the internal and external couplings 
are sufficiently weak, the smallest time scale in Eq. (41) is 0(1). Therefore, the explicit 
Euler scheme with the time step At ~ 0.01, which we used for the sake of simplicity 
and efficiency, is sufficiently accurate and stable under the parameter condition. 
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